Superfluid to Bose-glass transition of hard core bosons in one-dimensional 

incommensurate optical lattice 
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We study superfluid to Anderson insulator transition of strongly repulsive Bose gas in a one 
dimensional incommensurate optical lattice. In the hard core limit, the Bose- Fermi mapping allows 
us to deal with the system exactly by using the exact numerical method. Based on the Aubry- Andre 
model, we exploit the phase transition of the hard core boson system from superfluid phase with all 
the single particle states being extended to the Bose glass phase with all the single particle states 
being Anderson localized as the strength of the incommensurate potential increasing relative to 
the amplitude of hopping. We evaluate the superfluid fraction, the one particle density matrices, 
momentum distributions, the natural orbitals and their occupations. All of these quantities show 
that there exists a phase transition from superfluid to insulator in the system. 
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I. INTRODUCTION 

Anderson localization [l| was predicted fifty years ago 
as the localization of the electronic wave function in a 
disordered potential. So far the phenomenon of Ander- 
son localization has been observed in variety of systems, 
such as electromagnetic waves SB, sound waves [J], and 
quantum matter waves [E B & HI • Recently lots of at- 
tentions have been payed on the cold atom systems in 
disordered potentials [3, [l(| because of the experimental 
realization of the Anderson localization of the quantum 
matter waves 0, Q. The good tunability in ultracold 
atom systems offers myriad opportunities for study- 
ing the disorder effects in a controllable way. Random 
potential can be introduced by baser beams generating 
speckle patterns in the optical lattice [f| [l2j • By this way, 
Billy et al. observed exponential localization of a Bose- 
Einstein condensate (BEC) of 87 Rb atoms released into a 
one-dimensional waveguide 5]. Also one can achieve the 
random potential by loading a mixture of two kinds of 
atoms in a optical lattice with one heavy and one light, 
and the effect of the heavy ones is to produce an effective 
random potential for the lighter atoms when the heavy 
ones localize in the lattice randomly [1 31 ] . Another way to 
produce the random potential is to superimpose two one- 
dimensional (ID) optical lattices with incommensurate 
wave lengths to generate quasi-periodic potential [6|, . 
Particularly, Roati et al. @ observed localization of a 
noninteracting BEC of 39 K atoms in a ID system with 
two optical lattice potentials. 

On the other hand, the interactions between bosons 
in the cold atom systems can be controllably tuned by 
Feshbach resonances. It is thus feasible to experimen- 
tally study [HI the interplay between disorder and inter- 
actions in bosonic systems under controlled conditions, 
which has been a longstanding problem subjected to in- 
tensive studies. Theoretically, it has been predicted that 
there is a phase transition from a superfluid phase with 



extended sin gle part icle states to the Bose-glass (BG) 
phase [ll| Il7l.ll8llll| with Anderson localized single par- 
ticle states. The study of the interplay of disorder and in- 
teraction in strongly interac ting ultracold atomic systems 
also gets lots of attentionspaj2jJ^ Despite 
the intensive studies [H E^lllSlMllallE S3, Hi, 
the properties of the disordered interacting bosons re- 
main on debate. The theory of disordered interacting 
bosons is difficult and exact solutions are rarely known 
apart from numerical or approximate results. An excep- 
tion is in the one-dimensional (ID) limit with infinitely 
strong repulsive interactions, where the disordered boson 
problem can be exactly solvable via a Bose-Fermi map- 
ping [H. In the present work, we study the interact- 
ing bosons in the incommensurate optical lattices in the 
limit with infinitely repulsive interaction. In this limit, 
the Bose gas is known as Tonks- Girardeau gas or hard 
core boso nic g as [30| . which has been experimentally re- 
alized [3l|, 1.321 ]. The Bose-Fermi mapping or the Jordan- 
Winger transformation for hard-core bosons in ID op- 
tical lattice establishes a connection to non-interacting 
disordered fermions, which greatly simplify the compu- 
tation of physical quantities like the momentum distri- 
bution. Following an exact numerical approach proposed 
by Rigol and Muramatsu [HI], we calculate the super- 
fluid fractions, one particle density matrices, momentum 
distributions, natural orbitals and their occupations to 
exploit the phase transition of the system from super- 
fluid to BG phase. 



II. MODEL 

We consider N hard-core ultra-cold bosons in a ID 
incommensurate optical lattice. Under the single band 
tight binding approximation, the system can be described 
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by the following hard core bosons Hamiltonian: 
H = -tJ2(b\b %+1 + H.c.) + cos(a2ni + 5)n\ (1) 

i i 

where b\(bi) is the creation (annihilation) operator of the 
boson and they satisfy the hard core constraints [33j ] . 
i.e., the on-site anti-commutation ({bi,b\} — 1) and 
[bi,bj] = for i ^= j; n h i is the bosonic particle number 
operator; t is the amplitude of hopping and we will set 
it to be the unit of the energy (t — 1); V is the strength 
of incommensurate potential with a being an irrational 
number characterizing the degree of the incommensura- 
bility and S an arbitrary phase (in our calculation it is 
chosen to be zero for convenience, without loss of gener- 
ality). 

In order to obtain the ground state properties of the 
hard core bosons, we use the Jordan- Wigner transforma- 
tion 



particle wave function at i site with Vi the on site di- 
agonal potential and E n is the n-th single particle eigen 
energy with the eigen state given by |n) = J2i u i,nf}\fy- 
Obviously the solution to Eq. ([6j) is dependent on the po- 
tential Vi. For periodic Vi Bloch's theorem tells us that 
all the single particle states are extended band states. 
For a random Vi the system is the Anderson model and 
Eq.© produces [35| localized states in ID. For the in- 
commensurate case Vi = Vcos(27rm) with irrational a, 
Eq. ^ is the well-known Aubry- Andre model [36| . Aubry 
and Andre showed that when V < 2 all the single particle 
states are extended and when V > 2 all the single par- 
ticle states are the localized states. Mobility edges also 
have been found for the extended Aubry-Andre model 

MM. 



III. TRANSITION FROM SUPERFLUID TO 
BOSE-GLASS PHASE 



6 3 = f* n e ~ mslUi n e+ ™ flu fi > ( 2 ) 

0=1 13=1 

which maps the hard core bosons Hamiltonian into a non- 
interacting spinless fermions Hamiltonian 



H=- J^ififi+i + H.c.) + Vj2 C0S ( a2 ™)' 



(3) 



where fj(fi) is the creation (annihilation) operator of the 
spinless fermion and n{ is the particle number operator. 
The ground state wave function of the N spinless free 
fermionic system can be obtained by diagonalizing Eq. © 
and be written as: 



TV 



(4) 



n— 1 i—l 



where L is the number of the lattice sites, N is the num- 
ber of fermions (same as bosons), and the matrix P is 
given by the lowest N eigenf unctions of the Hamiltonian, 



P = 



( Pll Pl2 

P21 P22 



Pin \ 

P2N 



(5) 



V Pli Pl2 ■■■ Pln J 



with all the rows being ordered starting from the one 
with the lowest energy. 

The solution of Pj n reduces to a single particle prob- 
lem. The motion of a single particle in a ID lattice can be 
described by the nearest neighbor tight binding equation 



Ui+l,n + u i-l,n + ViUi 



(6) 



where the amplitude of the nearest neighbor hopping has 
been set to be unity, u^ n — Pi„ is the amplitude of the 



To see whether a bosonic system is in a superfluid 
phase, we can calculate the superfluid fraction, which 
reflects the response of a superfluid to the imposed phase 
gradient. Generally a nonzero superfluid fraction 
characterizes the system being in the superfluid phase. 
The superfluid fraction can be calculated with the aid 
of Peierls phase factors introduced in the Hamiltonian 
Eq.([T]) by means of the replacement bf&i+i — > b\bi + ie lv . 
In terms of free energy Fn(lp) for very small ip, the su- 
perfluid fraction is then determined as [3^, H(| 



Js tl f> N 



(7) 



For hard core bosons in one dimension, the above expres- 
sion can be represented as [4l| 



/f = 



L N 

i— 1 n— 1 

L N L 

~~k I y^[ M t,rt M »+l,™ 

n=N+lm=l i=l 

— u i+l, u i,m) I 3 

\N+1, 



(8) 



with UL+i, n = (— 1) u>i,n and E n is the eigenenergy 
of the n-th row in P matrix. The results of numerical 
calculations for the superfluid fractions for the hard core 
bosons in incommensurate lattices are shown in Fig[T]for 
three different lattice sizes and fillings. Here we con- 
sider odd numbers of bosons and use periodic boundary 
conditions [42j . When the strength of the incommensu- 
rate potential V is small, the system is in the superfluid 
phase with nonzero superfluid fraction. As V increases, 
the incommensurate potential acting as a pseudo-random 
potential makes the bosons difficult to hop and tend to 
be localized. Consequently, the superfluid fraction de- 
creases. When V gets around 2, the superfluid fraction 
approaches zero, and the system is in the Bose glass phase 
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FIG. 1: (Color online)The superfluid fraction versus V (the 
strength of the incommensurate potential) for three different 
lattice sizes and fillings with a — (y/E— l)/2. Insert: enlarge- 
ment of picture around the critical point. 
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FIG. 2: (Color online) Log-Log plots of the mean one particle 
density matrices for different V of the systems with 201 bosons 
on 1000 lattice sites and a = (v5 — l)/2. The red straight 
line shows a power-law fit — 0.195/\/i) for the mean one 
particle density matrix with V = 1.1. Insert: The log-line 
plot of the mean one particle density matrix for V = 2.1, 
and the red straight line shows an exponential-law fit (\np i = 
— 5.697 — 0.055i) for the mean one particle density matrix. 



(insulating phase) with all the bosons being localized [T^] . 
We note that with different irrational a and fillings the 
curves of superfluid fraction exhibit similar behavior. 

The one particle Green's function for the hard core 
bosons can be written in the form 



Gij — 



G 



HCB 



bibles) 



(9) 



\* G hcb) 



is the ground state of hard core 
bosons, and = (// nj»=\ e~ inf ^ \ = 



After a straightforward evalua- 



tion, the state \^ A ) can be written as 



JV+l L 



with P* = -Pi„ for I < i - 1, P£ = Pn for I > i with 
n < N, P A N+1 = 1 and P* +1 = (I ^ i), i.e., P A is 
formed by changing the sign of the elements P/„ for I < i 
because of the exponential term of the Jordan- Wigner 
transformation, and adding a column to P with element 
PiN+i = 1 and all the others equal to zero for the further 
creation of a particle at site i. In the same way we can get 
P B for the state \ 1 & B ) which has the same form as \^f A ) 
with the replacement of i by j. The Green's function is a 
determinant dependent on the L x (N + 1) matrices P A 
and P B (a! 



G, 



(V A \y B ) =det 



P A ) 1 P h 



(10) 



It follows that the one particle density matrix can be 
determined by the expression 

Pij = {b\b j ) = G ij +8 ij {l-2G ii ). (11) 
Then we define the mean one particle density matrix as 



Pi - T / ] Pj,j+i-l 

.7 = 1 



(12) 



to reduce the drastic oscillations caused by the incom- 
mensurate potential. The mean one particle density ma- 
trices for different V are shown in Figj2j The mean one 
particle density matrices for V < 1.1 are almost the same 
as the one with V = 1.1, having the power-law decay with 
exponent of —1/2, which has the same exponent of the 
hard core bosons in the lattice without the incommen- 
surate potential [33(. As V increases further, the mean 
one particle density matrix still has power-law decay, but 
the exponent is smaller than —1/2. Correspondingly the 
superfluid fraction decreases fast (see FigQ]). When V 
exceeds the critical point, the mean one particle density 
matrix has the exponential-law decay which is the char- 
acter of the Bose glass phase. 

The momentum distribution is defined by the Fourier 
transform with respect to i—j of the one particle density 
matrix with the form 



n(k) 



L 

2J = 1 



-ik(i-j) 



Pij, 



(13) 



where k denotes the momentum. In Fig. [3] we show the 
momentum distributions for systems with three differ- 
ent fillings in 1000 lattice sites. Odd numbers of par- 
ticles are taken for periodic boundary conditions. The 
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FIG. 3: Momentum distributions for systems with 1000 lattice 
sites, a = (VE - l)/2, V=l.l, and fillings of 201 (a), 501 
(b), 801 (c) hard core bosons, (d): The value of the tail of 
momentum distributions relates to the number of bosons (N) 
for systems with 1000 lattice sites, a = (y/E—l)/2 and V=l.l. 



peak structure in the momentum distributions reflects 
the bosonic nature of the particles, which is in contrast 
with the structure of the momentum distributions of the 
equivalent noninteracting fermions. Also because of the 
existence of the additional lattice to produce the incom- 
mensurate potential, it is possible that there exist other 
peaks besides at k = in the momentum distributions. 
For a — (V5 — l)/2, the secondary peaks are found to 
appear at k = ±(a — l)2n only for low and high fillings. 
On increasing the number of the bosons the population 
of the high momenta states is always increasing, which 
leads to the increase in the tails of the momentum dis- 
tribution as shown in Fig[3ji. 

Then we consider the properties of the momentum dis- 
tributions as the strength of the incommensurate poten- 
tial changing. In Fig[5] we show the momentum distri- 
butions for systems with three different V. Low fillings 
are required for the existence of secondary peaks in the 
momentum distributions. On increasing V , the value of 
n,k=o decreases which means that the coherence among 
the hard core bosons decreases and the superfluid frac- 
tion decreases. Also the peak at k = and the sec- 
ondary peaks all become widespread. The value of the 
secondary peaks (FigfJJl) first increases, and it starts to 
decrease when V is bigger than around 1.5. Finally it 
reaches an almost fixed number as the system going into 
the BG phase. The population of high momenta states is 
always increasing accompanying the decrease of the peak 
at k = 0. 

Now we consider the properties of the secondary peaks. 
Their positions are only decided by the value of a and 
are irrelevant with the system size, V and filling. But 
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FIG. 4: Momentum distributions for systems with 1000 lattice 
sites, a = (\/5-l)/2, 201 hard core bosons, and V = 1.1 (a); 
V = 1.7 (b); V = 2.3 (c). (d): The value of secondary peaks 
of momentum distributions relates to V for systems with 1000 
lattice sites, 201 bosons and a = (v5— 1)/2. (e): The position 
of the secondary peak (the k > one) in the momentum 
distributions relates to a' for systems with 1000 lattice sites, 
201 hard core bosons, and V = 1.1. The irrational numbers 
beside the points are the corresponding a. 



the existence of the secondary peaks are related to the 
filling, and the peaks only exist for low and high fillings. 
For the incommensurate lattice with the form cos(27rm), 
actually we can restrict a in the range of [0, 0.5] because 
cos[27r(a— = cos(27rm) and cos(— 27rcn) = cos(27rm). 
For example, a = (yo± 1) /2 are the same for the system. 
So any a € (— oo, oo) has an equivalent number in the 
range [0, 0.5], and we denote it by a'. The positions of the 
secondary peaks are decided by the value of a'. In FigfJJ: 
we show the position of the secondary peak (the k > 
one) as a function of a! . From the data we can see that 
the positions of the secondary peaks in the momentum 
distribution are decided by a and at k = ±2a'7r. 

The natural orbitals are defined as the eigenfunc- 
tions of the one particle density matrix [43[, 
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FIG. 5: Occupations of the natural orbitals for systems with 
1000 lattice sites, 201 bosons, a = (\/E - l)/2 and V = 
(a); V = 1.5 (b); V = 2.3 (c). (d): The amplitude of the 
discontinuation (Z) defined as Z = Ajv — Ajv+i relates to 
V for systems with 1000 lattice sites, 201 bosons and a — 
(V5-l)/2. 
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FIG. 6: Profiles of the two lowest natural orbitals for systems 
with 1000 lattice sites, 201 bosons, a = (\/E- l)/2 and V = 
(a)(b); V = 1.5 (c)(d); V = 2.3 (e)(f). 



and can be understood as being effective one particle 
states with occupations A^. For noninteracting bosons, 
all the particles occupy in the lowest natural orbital 
and bosons are in the BEC phase at zero temperature 
(only the quasi condensation exists for the ID hard core 
bosons). The occupations of the natural orbitals for sys- 
tems with three different V are shown in FigOO The 
occupations are plotted versus the orbital numbers rj, 
and ordered starting from the highest occupied one. For 
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FIG. 7: The condensate fraction f c versus V for systems with 
1000 lattice sites, 201 bosons, and a = (V5- l)/2. 



V < 2, the occupation distribution exhibits sharp single- 
peak structure. The peak appearing in the lowest orbital 
is the feature of the boson which is against the step func- 
tion of the fermions. With the increase in the strength 
of incommensurate potential, the occupation of the low- 
est natural orbital (Ai) decreases. When V > 2, no an 
obvious peak appears in the lowest natural orbital. We 
also find that a discontinuation at rj = N emerges when 

V > 2. To characterize such a discontinuation, we define 
Z = Xjst ~ Xn+i which indicates the occupation differ- 
ence between the N-th. and (N + l)-th natural orbital 
for a boson system with N particles. The amplitude of 
the discontinuation Z versus the strength of the incom- 
mensurate potential (V) is plotted in Fig[5ji. There is 
an obvious change around V = 2. For V < 2 there is 
no discontinuation in the occupation number. However, 
for V > 2 a nonzero Z appears and the amplitude of the 
discontinuation increases with the increase in V. 

The effect of incommensurate potential on the natural 
orbital is shown in Figj6l where profiles of the two lowest 
natural orbitals for the same system with three different 

V are plotted. When V = 0, there is only the periodic 
optical lattice, and the natural orbitals are plane waves. 
As V increases but is smaller than 2, the natural orbitals 
still spread over all the lattice corresponding to extended 
states, but there are a lot oscillations in the waves in- 
duced by the existence of the incommensurate potential 
which acts like random potential on sites because of the 
irrational a. When V > 2 the states do not spread over 
all the lattice any more and are localized. Correspond- 
ingly the system is in the BG phase. 

Finally we consider the influence of incommensurate 
potential on the condensate fraction, which is defined as 
fc = Ai/iV to indicate the ratio of occupation of lowest 
natural orbital. In Figl7]we show the condensate fraction 
as function of V. For small V the condensate fraction f c 
decreases slowly with the increase m V. As V increases 
further to approach V = 2, the condensate fraction de- 
creases rapidly. For V > 2 there is almost no conden- 
sation. The change of the condensate fraction also gives 
signature of the superfluid to insulator transition in the 
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incommensurate optical lattice system. 



IV. SUMMARY 

In summary, we have studied the properties of hard 
core bosons in an incommensurate optical lattice. Us- 
ing the Bose-Fermi mapping and the exact numerical 
method proposed by Rigol and Muramatsu [33[, we ex- 
ploit the phase transition from superfluid to the localized 
BG phase as the strength of the incommensurate poten- 
tial increases from weak to strong. We calculate the su- 
perfluid fraction, one particle density matrices, momen- 
tum distributions, the natural orbitals and their occu- 
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